{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Boson Sampling with MPS\n", "\n", "In this notebook, we explain how to use the MPS (Matrix Product State) backend to simulate a linear circuit. MPS simulation is based on a type of tensor network simulation, which gives an approximation of the output states [[1]](#References) [[2]](#References). It does the computation on each component of the circuits one-by-one, and not on the whole unitary. The states are represented by tensors, which are then updated at each component. These tensors can be seen as a big set of matrices, and the approximation is done by chosing the dimension of these matrices, called the *bond dimension*. For this example, we simulate a simple boson sampling problem, with 6 modes and 3 photons [[3]](#References)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import perceval as pcvl\n", "import matplotlib.pyplot as plt" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Definition of the circuit\n", "\n", "Just like in the Boson Sampling notebook, we generate a Haar-random unitary and its decomposition in a circuit :" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6 modes triangular Boson Sampler :\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Φ=1.265255\n", "\n", "\n", "Φ=2.858772\n", "\n", "\n", "Φ=3.402583\n", "\n", "\n", "Φ=1.851956\n", "\n", "\n", "Φ=2.413235\n", "\n", "\n", "Φ=0.39409\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=2.758766\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=3.479159\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.947194\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=2.84937\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.169331\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=2.251856\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.194179\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.333985\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.036737\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=2.642193\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=2.172429\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=2.095071\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.661671\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.738002\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.273074\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=3.543207\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.449777\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=3.302943\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.02537\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.249399\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=6.238378\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=2.976814\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.71687\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=3.793155\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.998835\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.329839\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.382743\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.389795\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.488344\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=3.753082\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "0\n", "1\n", "2\n", "3\n", "4\n", "5\n", "0\n", "1\n", "2\n", "3\n", "4\n", "5\n", "" ], "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = 6\n", "unitary = pcvl.Matrix.random_unitary(m)\n", "mzi = (pcvl.BS() // (0, pcvl.PS(phi=pcvl.Parameter(\"φ_a\")))\n", " // pcvl.BS() // (1, pcvl.PS(phi=pcvl.Parameter(\"φ_b\"))))\n", "linear_circuit = pcvl.Circuit.decomposition(unitary, mzi,\n", " phase_shifter_fn=pcvl.PS,\n", " shape=\"triangle\")\n", "\n", "print(m, \" modes triangular Boson Sampler :\")\n", "pcvl.pdisplay(linear_circuit, compact = True)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## MPS simulation\n", "\n", "Let us now define the MPS simulation, using the MPS backend in Perceval." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "mps = pcvl.MPSBackend()\n", "mps.set_circuit(linear_circuit)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We now choose the size of the matrices (the bond dimension) for our simulation." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "chi = 8\n", "mps.set_cutoff(chi)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "And finally, we get the output probability distribution from a given input state with 3 photons." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
state probability
|0,0,1,0,2,0> 0.089692
|1,1,1,0,0,0> 0.074157
|0,1,0,0,0,2> 0.066353
|0,0,1,0,1,1> 0.057677
|0,1,1,0,0,1> 0.057048
|0,0,0,0,1,2> 0.052347
|0,1,1,0,1,0> 0.048586
|0,1,2,0,0,0> 0.039534
|1,1,0,0,0,1> 0.034403
|2,0,1,0,0,0> 0.034118
|1,0,2,0,0,0> 0.029711
|0,0,0,1,2,0> 0.029622
|0,0,0,0,3,0> 0.027699
|0,0,2,0,1,0> 0.026586
|0,0,0,1,1,1> 0.022117
|2,0,0,0,0,1> 0.021964
|0,1,0,1,0,1> 0.020708
|1,1,0,1,0,0> 0.020591
|1,0,0,0,0,2> 0.018927
|0,0,2,0,0,1> 0.01776
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n = 3\n", "input_state = pcvl.BasicState([1]*n + [0]*(m-n))\n", "mps.set_input_state(input_state)\n", "probs = mps.prob_distribution()\n", "pcvl.pdisplay(probs, max_v=20)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Certification of the MPS method :\n", "\n", "As we make an approximation by chosing the bond dimension, we have to check when does this approximation becomes good enough. Unfortunately, there is no formula giving the minimal size for a given approximation error. What we can do though is to compute the *Total Variance Distance* (TVD) between an ideal simulation of Boson Sampling, and an approximated one. To compute the ideal one, we can for instance use the *SLOS* backend on Perceval :" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
state probability
|0,0,1,0,2,0> 0.089692
|1,1,1,0,0,0> 0.074157
|0,1,0,0,0,2> 0.066353
|0,0,1,0,1,1> 0.057677
|0,1,1,0,0,1> 0.057048
|0,0,0,0,1,2> 0.052347
|0,1,1,0,1,0> 0.048586
|0,1,2,0,0,0> 0.039534
|1,1,0,0,0,1> 0.034403
|2,0,1,0,0,0> 0.034118
|1,0,2,0,0,0> 0.029711
|0,0,0,1,2,0> 0.029622
|0,0,0,0,3,0> 0.027699
|0,0,2,0,1,0> 0.026586
|0,0,0,1,1,1> 0.022117
|2,0,0,0,0,1> 0.021964
|0,1,0,1,0,1> 0.020708
|1,1,0,1,0,0> 0.020591
|1,0,0,0,0,2> 0.018927
|0,0,2,0,0,1> 0.01776
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "slos = pcvl.SLOSBackend()\n", "slos.set_circuit(pcvl.Unitary(unitary))\n", "slos.set_input_state(input_state)\n", "probs_slos = slos.prob_distribution()\n", "pcvl.pdisplay(probs_slos, max_v=20)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We also have to define the TVD function." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def tvd(probs1, probs2):\n", " tvd = 0\n", " for state, prob in probs1.items():\n", " tvd += abs(prob - probs2[state])\n", " return tvd\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now we compute the TVD between the two simulations for different bond dimensions, going from 1 to 10." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAEWCAYAAADRrhi8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4hElEQVR4nO3dd3wUdfrA8c+zqZCEmtBLEEITRYRTBA9RUNE7e8XTs2I/Pevp6c+fp573O3vDw3rYu57c2cVOUQGxULIECE0SQt8QIO35/TETbonZFLK7s5t93q9XXtmdmZ15ZnZnnvnOfOf7FVXFGGOMiSc+rwMwxhhjmsqSlzHGmLhjycsYY0zcseRljDEm7ljyMsYYE3cseRljjIk7niUvEVER6beHn/21iOSHO6ZGLHeAiMwXkYCIXNHIz+zxenpNRM4Rka+8jsM0jYhMFZE7vI4jUpp57CgUkfEhxu12XAmeVkT+LCJP7lnEDcaU665Tsvv+PRE5O0zzDrlOYZr/AhEZG675NUWDyctd2e0iUhr090g0gguKYbcfq6p+qaoDohmD63rgU1XNUtWHao8Ukc9E5IJwL9T9AdZs+23u9gj+Pp4WkWfr+NxQEdkpIh1E5FYRqXATb0BE/CLyiIh0DXe87rJjPvGJSDt32xUFbZMbgsaHPEiKyG9F5Bv3+9ggIi+ISI+g8akicq+IrHa/o0IReSAKq1Uv93tREbm/1vDj3OFT3fc1B9TSoPhvqDX9fBHZKiLrReQTEekT5dUJq/qOK6p6p6peAL9MNhGI4yhVfaah6RqTxMN5rKzrpEhV91bVz8Ix/6ZqbMnrGFXNDPq7PKJRxa7ewIJoL9T9AWaqaiawtzu4XdCwJ4ATRSSj1kfPAv6jqhvd96+oahbQATgB6ALMjVQCiwP3A5nAIKAtcCxQ0NCHRORk4EXgASAb5zvZCXwlIu3dyW4ERgAHAFnAWGBeWKPfc0uBU2sdfM8G/HVM2879jU0EbhGRCe4B81ngGpzt1geYDFQ1N7BIJYRE1OK3parW+wcUAuPrGJ4GbAaGBA3LAbYDndz3k3AOBhuBaUC3oGkV6Oe+/gy4IGjcOcBX7usv3Gm3AaXAaTgHgtVB0w9y57EZJ7kcGzRuKs6O9Q4QAL4G+tazvse689jsznOQO/wTnJ1zhxtH/1qf+2ut8Y8ErefFwBJ3npMBCfrcecAiYBPwAdC7ge8j151ncq3h+cDvg94nAT8Dx7nvbwWer/WZJOB74J4QyzoHmAE8AmwBFgPjgsa3BZ4C1gJrgDvceQ5yt0OVuy024xzgNgM+97NPAOuC5vUc8Mf65tuYbdbQ9q61fj8Bx9ezrXf9RoOGCbACuL7WcJ87v9vc9/+pWZ/G/AEPAquArcBc4NdB424FXsVJGAGc3+eIoPHDcBJjAHgFeBm4o57v9CvgfeA37rAOQBFwNzA11O8M+Ba4FjgZmN+EdZsKTAE+cmP8vI7v7DL3O1veyGPHFcAyYL0bd83vqi/OvrrBHfcCTgIOPp7dCCx0fz//BNLdcWPZ/bhSiHvsI2j/AVa6MZS6f4e4ce4T9NlOQBmQU8f2SALuceNb5q77rm1N0PEQ6Odury3u9K80dFwE/uR+n8+FWKdQ638O7nG39j4AXAhUAOXu8v5dxzZKwzmh+9n9ewBIC962OCc863D27XODlnO0G1MAZ5+/tsHfVSN+eLuCq2Pc08Bfg95fBrzvvj7M3dj7uyv1MPBFXQcG6kledR1Egr8QIAXnR/5nINVdbgAYELTjbMA5A07G+TG/HGJ9+rs/hsPd+V7vzju1rjjr+Pwvxrux/wdoB/QCSoAJ7rjj3PkPcmO7GZjZwPeRS93J6ybg46D3R7rLSqm989X63G3A1/Uc6CqBq9ztcRrOTtTBHf8W8BiQgbOzfgNcVM+OsBIY7r7Ox9lxBwWNG9aI+da7zerb3nWs35M4ieBcIK+O8XUlr4Hu8D51TP8XYJb7+mZ3nS4F9iFEAg367JlAR3edrsE5+NQcVG7FORk4GufA9zdgtjsuFSeZ1nxHJ+McZBpKXmfw3wPhpe72voM6khdOwh6NczAeB+zlxnM/cCiQ2cC6TcXZJ8fgHAse5Jf790c4SbQVjTt2fOpO3wunxBh8sD/c/VwOzkH+gVrHs5+Anu7nZ9RsKxqfvHZtm6BpHwX+HvT+StwDfB3b42KcE8GaGD4ldPJ6CWff9gHpwMENHBcrgb+7698qxDqFWv9zCJG8gr7HO2qND95GtwGzcfbZHGAmcHut2G7D+Z0ejfN7au+OX4t7wga0B/av7zel2vjkVXP2XPM3yR03HlgaNO0M3LN/nDPnu4LGZeLsVLl1bJRdX1ZdGzHEl1STvH6Ns6P7gsa/BNwatMGfDBp3NLA4xLr+D/Bq0HsfzlnA2LrirOPzvxjvxh78g3sVuMF9/R5wfq3llVFP6YvQyauXu317uO9fAB4MGn8rdSevi4ElIZZ1Ds4ZVHBJ8Rucy5GdcS6VtQoaNxHnnuAvvkN32HPA1TiXK/OBu9zl7yqVNWK+9W6z+rZ3HevXCuekZ6677QqAo0L97txhB7vD0+vbljhJ5jKcfWKnux3Pbmh/C5rXJmBo0HcXfGIyGNjuvh5Tx3c0k4aTVyugGKeUOxsnOdWVvDa7sSwCrgiaz0h325bgJLKphEhi7riXg95n4pTKewZt58OCxjfm2DEhaPylwPQQyz4e+C7ofSFwcdD7o3GPYTQveR2Ic7Ii7vs5wKkhYvqkVgxHEDp5PQs8jrtf15pPXcfFcoJ+myHWKdT6n0PzktdS4OigcUcChUFxbK+1zdYBI93XK4GLgDaN3Ucae8/reFVtF/T3hDv8U6C1iBwoIrnAfjhnzQDdcM4IAVDVUpwSUPdGLrOxugGrVLU6aNiKWsspCnpdhrMzhJpXcMzVOJdymhtzqOX3Bh4Ukc0ishnn0oPsyfJUdSXOWeaZIpKJs9P+ohJHHbq7yw1ljbq/LtcKnO3UG+cMam1Q/I/hnHWF8jnOj3iMG+tnOJdcDgG+dLd3Q/NtzDZr1PetqtvVuRE/HKfU8yrwmoh0qGcd1rv/67pP2LVmvKpWqepkVR2NUwr8K/C0iAyqa6Yicq2ILBKRLe56tcW5nxZqndLdexrdqPs7qpeqbse5lH4z0FFVZ4SYNFtV26vqIA2qpKSqs1X1VFXNwTmBHINTQghlVdBnS3G+t251jadxx47g6Wt+k4hIZxF5WUTWiMhW4Hl2344hP9scqvo1zvcyVkQG4pQAp4WYvFsdMYRyPc7v+xu3Zt95DYRSoqo7Gpgm7Ovv2u17q2PeG1S1Muh98L55Ek4iXSEin4vIQQ0trFlV5VW1CmeHn+j+/UdVA+7on3EONAC4lQk64pRkatsGtA5636UJYfwM9BSR4HXpFWI5jZlXcMyCU7xu7Ly04Ul2swrncljwiUErVZ3ZxPnUeAanVHQSzr2DufVN7G6zY4Av65msu7sdavTC2U6rcEoU2UGxt1HVmgoldW2Lz3EOdGPd11/hnPEf4r6nEfMN9zZzglXdCtyJc6myvlpz+TjX7k8JHuhuy5OA6XXMe7uqTsYpwQyuPV5Efo1zkDoV5zJKO5zLs1J72jqspe7vqDFqKl0838jp66Sq3wJvAkPqmaxnzQv35KoDzu9o12yCXjfm2NEz6HXNbxKc71Bx7j+1wbkcW3s7hvpsY4Xaz59xl3cW8Ho9SWRtHTHUvSDVIlWdpKrdcEomjzZQw7Axx6BQ67/bcVhEah+HG5r3bt8bTdi2qvqtqh6Hc5L6L5y8Uq9wPOf1Is69kN+5r2u8BJwrIvuJSBrOj+prVS2sYx7zcWrLtXa/mPNrjS/Guc5el5oznutFJEWcZw6Owblp3VSvAr8RkXEikoKzY+/EuQzTGPXFWZcpwI0isjeAiLQVkVMa+Ex93sD5wfwFZ0eqk4gkuyWAl3BOFO6rZ56dgCvcbXsKzr2md1V1LfAhcK+ItBERn4j0FZFD3M8VAz1EJLVmRqq6BOfSwZnA527CKMY56H/uTtPQfMO2zUTkf0TkV2619nSc+xSbcRJUjVQRSa/5w9lnrgVuFpEz3OFdcO6ftcG5D4SI/FFExopIK3d7n41T6/C7OkLJwrkfUAIki8gt7rwaY5b72Zrv6ESc+7uN8TnO/aGHGzk9ACJysIhMEpFO7vuBOBWdZtfzsaPdz6UCt+Pcs1sVYtrGHDuuE5H2ItIT53t7xR2ehXObY4uIdAeuq2P+l4lID7eEfVPQZxurBKjml/v68zi1eM+k/qser+J8Xz3c2qk3hJpQRE6R/z6CsQkngdRcZWrq8aZGqPX/Htjb3e7pOJdKgzW0vJdw9oscEckGbqERJ0bu/vc7EWmrqhU4lZaqG/pcY5PXv2X354pqLg3WFJe34RQP3wsa/jHOPaQ3cM40+gKnh5j//TjXaotxDrov1Bp/K/CMe6no1OARqlqOk6yOwrlk8yjOfbfFjVy34Hnl4/zwHnbndQzOYwLljZzFg8DJIrJJRH7xHFgdy3sL5+bqy+4ljp/c9dgjqroNZ3v34JfbEOA0ESnFOaufhnMpZriq1nd29DWQh7M9/gqcrKob3HG/x6kwUFNz6XX+ezntE5zKEEUisj5ofp/jXD5YFfRe2L0aecj5hnmbKU5tq/U4Z4iH49TAKw2aZgFOwq35O1dVX8E5u74KZxsuxLmHNDpo25QB9+Jc7luPc//rJFVdVkccH+DU/vPjXGrZwe6XdkKvgPPbPBHnfsVGnBPJNxv5WVXV6frfRykaazNOsvrR/T29j3O74K56PvMi8L9ujMNx9rNQcTXm2PE2zr3K+TiXP59yh/8Fp6LHFnd4XdviRZwTpGU492ma9EC3qpbh7Asz3GPSSHf4KpzfsVL/1YwncL7z793p6/u+fgV87W7nacCVQb+hWwlxXGxAneuvqn6cChUf49T8rP2c5lPAYHd5/6pjvnfg3Ov7AfjRXbfGbtuzgEJ3n74YpzBUr5qbi8YYExHiPPi8WlVv9jqWSBORp4GfE2FdvdayH2IzxpgoEafS2ok4z92ZCLOGeY0xpplE5HacS9h3q+pyr+NJBFG7bOgWp3+L06rCL2olubWlHuS/D6+do6qx0pyOMcaYGBLNktdUYEI944/CqRiQh9MUyT+iEJMxxpg4FLV7Xqr6hXtNOJTjgGfdhy1ni9Pid1e36nRI2dnZmptb32yNMcbUNnfu3PXuQ+ZxKZYqbHRn9+rBq91hv0heInIhTumMXr16MWfOnKgEaIwxLYWINNgSSyyLywobqvq4qo5Q1RE5OXF74mCMMWYPxVLyWsPuzZb0YM+aeDLGGNPCxVLymgb8XhwjgS0N3e8yxhiTmKJ2z0tEXsJpkDVbRFbjNBWTAqCqU4B3carJF+BUlT83WrEZY4yJL9GsbTixgfGK0/6bMcYYU69YumxojDHGNIolL2OMMXHHkpfHlpaU8tqcVZTurGx4YmOMMYAlL8898PESrnv9B0beOZ1bpy1gWUlpwx8yxpgEF0stbCSk/KKtDO3Zjr7ZGbzw9QqmzixkTP8czhnVm7H9O+HzNaYneGOMSSxW8vJQeWU1y0q2cXC/jtx32n7MvGEcVx/en/yirZw3dQ6H3vsZT365jC3bK7wO1RhjYoolLw8tW19KZbXSv3MWADlZaVwxLo+v/nQYD08cRk5mGne8s4iRd07nprd+xF8c8DhiY4yJDXbZ0EP5RU4yGtAla7fhKUk+jhnajWOGduOnNVt4ZmYhr81dzQtfr2RU346cPSqX8YM6k2SXFI0xCcpKXh7yFwdI9gl7ZWeGnGZI97bcfcpQZt84jusnDKBw/TYuem4uY+76lCmfL2XTtvIoRmyMMbHBkpeH8otK6ZOdQWpyw19Dh4xULh3bjy+uP5QpZ+5Pzw6t+L/3FjPyb9P50+s/sPDnrVGI2BhjYoNdNvSQvzjAPj3aNukzyUk+JgzpyoQhXVlctJVnZq7gre9W88qcVRyQ24GzR+VyxN6dSUmy8xJjTMtlRziPlJVXsnJjGQM7ZzU8cQgDu7Thbyfuw9c3juemowexdut2LntxHr/++6c88skS1pfuDGPExhgTO6zk5ZElxc7DyP277HnyqtG2dQqTxuzFeQf34dPF63hmViH3fOjnoekF/HZoV84Zlcu+Pdo1eznGGBMrLHl5JN+t9j6gGSWv2pJ8wvjBnRk/uDMF60p5dlYhb8xdzZvz1jCsVzvOGZXLUUO6NuoemzHGxDI7innEXxQgPcVHzw6tIzL/fp0yue24Icz68zj+95jBbC6r4MqX5zP6759w/0d+1m3dEZHlGmNMNFjJyyP5xQHyOmVF/FmtNukpnDu6D2cflMsXS0p4ZmYhD05fwqOfFXDUkK6cPSqX/Xu1Q8SeGTPGxA9LXh7xFwc4uF9O1Jbn8wljB3Ri7IBOLF+/jWdnFfL6nNVM+/5n9unelrNH5XLM0K6kJSdFLSZjjNlTdtnQA5vLyineupMBXUI/nBxJfbIz+N9j9mb2n8dx+3F7s72iimtf+57TH5/NjooqT2IyxpimsOTlAX9NTcMwVtbYExlpyZx1UC4fXTWGe08ZyncrN3PTWz+hqp7GZYwxDbHLhh7YVdMwDNXkw0FEOGl4D1ZtKuOBj5cwuFsbzj+4j9dhGWNMSFby8kB+0Vay0pPp0ibd61B2c8VheRy5d2f++s5Cvlqy3utwjDEmJEteHvAXlTKgc1bM1fDz+YR7T92PvE5ZXPbiPFZs2OZ1SMYYUydLXlGmquQXB8LSskYkZKYl88TvRyACk56dQ+nOSq9DMsaYX7DkFWXrAjvZsr0irC1rhFuvjq2ZfMb+LC3ZxjWvzqe62ipwGGNiiyWvKKvpgNLrmoYNGd0vm5uOHsQHC4p56JMlXodjjDG7sdqGUeaPsZqG9Tl3dC4L127lgY+XMLBLGyYM6eJ1SMYYA1jJK+ryiwLkZKXRISPV61AaJCLccfwQ9uvZjqtfnc/iIuvw0hgTGyx5RZm/OBDT97tqS09J4rGzhpOZlsykZ+ewaVu51yEZY4wlr2iqrlb8xaUxf7+rts5t0ply1nCKt+zk8pfmUVlV7XVIxpgEZ8krilZv2s72iirP2jRsjv17teeOE4Ywo2ADd7672OtwjDEJzipsRFFNs1DxVvKqceqInixau5WnZyxnUNcsThnR0+uQjDEJykpeUVRT0zAvTpMXwE1HD2JU347c9NZPfLdyk9fhGGMSVFSTl4hMEJF8ESkQkRvqGN9LRD4Vke9E5AcROTqa8UVaflGAHu1bkZkWvwXe5CQfk8/Yn85t07joubkUW4/MxhgPRC15iUgSMBk4ChgMTBSRwbUmuxl4VVWHAacDj0YrvmiIt5qGobTPSOWJ34+gdGclFz031/oAM8ZEXTRLXgcABaq6TFXLgZeB42pNo0Ab93Vb4OcoxhdRFVXVLC0pjdk2DZtqYJc23HfqUOav2szN/7I+wIwx0RXN5NUdWBX0frU7LNitwJkishp4F/hDXTMSkQtFZI6IzCkpKYlErGFXuH4bFVXaIkpeNSYM6cqV4/J4fe5qps4s9DocY0wCibUKGxOBqaraAzgaeE5EfhGjqj6uqiNUdUROTk7Ug9wTi+OkTcOmunJcHkcM7swd7yxiRoH1AWaMiY5oJq81QHDd6h7usGDnA68CqOosIB3Ijkp0EeYvDpDkE/bKyfA6lLDy+YT7TtuPvjkZXPbiPFZuKPM6JGNMAohm8voWyBORPiKSilMhY1qtaVYC4wBEZBBO8oqP64INyC8KkNuxNekpSV6HEnY1fYCpOn2AbbM+wIwxERa15KWqlcDlwAfAIpxahQtE5DYROdad7Bpgkoh8D7wEnKMtpCaAvzgQFy3J76neHTOYfMb+LFkX4GrrA8wYE2FRveelqu+qan9V7auqf3WH3aKq09zXC1V1tKoOVdX9VPXDaMYXKdvLq1ixsazF3e+q7eC8bP7s9gH28CcFXodjjGnBYq3CRotUsK4UVVpUTcNQzj+4Dyfu3537P/bzwYIir8MxxrRQlryiID+OOqBsLhHhzhP2YWiPtlz9yvxdPUcbY0w4WfKKAn9xgNRkH707tqyahqE4fYCNoLXbB9jmMusDzBgTXpa8oiC/KEBep0ySfOJ1KFHTpW06U84cTtGWHVz+4nfWB5gxJqwseUVBS2nTsKmG93b6APuqYD1/e8/6ADPGhE/8Nm8eJ7Zsr2Dtlh0tpk3Dpjp1RE8W/ryVp75azuCubThpeA+vQzLGtABW8oqwJTWVNRKw5FXjpt84fYDd+NaPzF+12etwjDEtgCWvCNvVe3KClrwAUmr6AGuTxkXPzWGd9QFmjGkmS14R5i8KkJmWTLe26V6H4qmaPsACOyq56Pm57Ky0PsCMMXvOkleE5RcH6N85E5HEqWkYysAubbj3lKF8t3IzN79lfYAZY/acJa8IUlXyi1p2m4ZNddQ+XbliXB6vzV3NM9YHmDFmD1nyiqCS0p1sKqto8W0aNtUfx+Vx+ODO3P7OImZaH2DGmD1gySuC/EWlQGLXNKyLzyfcd+pQ9srO4NIX57Fqo/UBZoxpGkteEWQ1DUPLSk/hid+PoLparQ8wY0yTWfKKIH9RgI4ZqWRnpnkdSkzKzc7gkTP2x18c4NrXvrc+wIwxjWbJK4KcmoZW6qrPmP45/PnoQbz3UxGPfGp9gBljGseSV4RUVytLWnjvyeFy/sF9OHFYd+77yM+H1geYMaYRLHlFyJrN29lWXmXJqxFEhDtPdPoAu+qV+Swu2up1SMaYGGfJK0L8NZU17LJho9T0AZaRlszpj89m3spNXodkjIlhlrwiZFdNw86ZHkcSP7q0Tee1iw+ibasUznhiNp8sLvY6JGNMjLLkFSH+ogDd27UiKz3F61DiSu+OGbx+8Sj6dcpk0rNzeW3OKq9DMsbEIEteEZJfXGqlrj2Uk5XGyxcexEF7deS613/g0c8KrB1EY8xuLHlFQGVVNUvXldrDyc2QmZbM0+f8imOHduOu9/P5y78X2nNgxphdrCflCCjcUEZ5VbU1C9VMqck+HjhtP7Iz03h6xnLWl+7k3lOHkpac5HVoxhiPWfKKAKtpGD4+n/A/vx1EpzZp/N97i9lUVs6UM4fbvURjEpxdNoyAxUUBfAL9Otk9r3AQES4+pC/3njKU2cs2cvrjsykJ7PQ6LGOMhyx5RYC/KEBuxwzSU+zyVjidNLwHT549gmUl2zjpHzMpXL/N65CMMR6x5BUBfmvTMGIOHdCJFycdSGBHBSdPmcmPq7d4HZIxxgOWvMJsR0UVhRu2WU3DCBrWqz2vXzKKtOQkTn98Fl8uKfE6JGNMlFnyCrOCdaVUq3VAGWl9czJ589JR9OzQmvOmfsvb89d4HZIxJooseYVZTU3DAV2sskakdW6TzisXHcSwXu258uX5PPXVcq9DMsZEiSWvMMsvDpCa5CO3Y4bXoSSEtq1SePa8A5iwdxdu/89C/u+9xdYahzEJwJJXmPmLAvTtlElykm3aaElPSWLy7/bndwf2YsrnS7nmte+pqKr2OixjTARF9QgrIhNEJF9ECkTkhhDTnCoiC0VkgYi8GM34wsFfXMoAa9Mw6pJ8wh3HD+Hqw/vz5rw1THp2DmXllV6HZYyJkKglLxFJAiYDRwGDgYkiMrjWNHnAjcBoVd0b+GO04guHwI4K1mzebjUNPSIiXDEujztP2Icv/CWc8cTXbNxW7nVYxpgIiGbJ6wCgQFWXqWo58DJwXK1pJgGTVXUTgKqui2J8zeYvLgWspqHXzjiwF/84cziL1m7l5CkzWb2pzOuQjDFhFs3k1R0I7pxptTssWH+gv4jMEJHZIjKhrhmJyIUiMkdE5pSUxM4zPtamYew4cu8uPH/BgawP7OTER2eyaO1Wr0MyxoRRrNUqSAbygLHAROAJEWlXeyJVfVxVR6jqiJycnOhGWI/8ogAZqUl0b9fK61AM8KvcDrx28Sh8Ipz62CxmL9vgdUjGmDCJZvJaA/QMet/DHRZsNTBNVStUdTngx0lmccFfHCCvcxY+n3gdinEN6JLFG5eOolNWGr9/+hve/2mt1yEZY8IgmsnrWyBPRPqISCpwOjCt1jT/wil1ISLZOJcRl0UxxmbxFwfsflcM6t6uFa9fPIq9u7Xhkhfm8fzsFV6HZIxppqglL1WtBC4HPgAWAa+q6gIRuU1EjnUn+wDYICILgU+B61Q1Lq71rC/dyfrScqtpGKPaZ6Ty4gUjOWxAJ27+10/c95HfHmY2Jo5FtTNKVX0XeLfWsFuCXitwtfsXV/xFbrNQVvKKWa1Sk3jsrOHc+OaPPDR9CSWBndx+3N72QLkxcch6Ug6T/JqahtamYUxLTvJx18n70qlNGpM/XcqG0p08NHGY9b1mTJyxU84w8RcHaN86hZzMNK9DMQ0QEa47ciC3HjOYjxYVc9ZTX7OlrMLrsIwxTWDJK0zyi5wOKEWspmG8OGd0Hx6eOIz5qzZz6mOzWLtlu9chGWMayZJXGKiq06ahVdaIO7/dtxvPnHsAazZv56RHZ1KwLuB1SMaYRrDkFQY/b9lB6c5Ka1kjTo3ql83LF46kvEo5ecos5q7Y5HVIxpgGWPIKg5qahgOt5BW3hnRvy5uXjKJdqxR+9+Rspi8q9jokY0w9LHmFQU1NwzwrecW1Xh1b8/olo8jrlMWFz81lZsF6r0MyxoRgySsM/EUBurZNp22rFK9DMc2UnZnGSxeOpGvbdO58b5E9yGxMjLLkFQb5xQG739WCZKYl88fx/flpzVbe/6nI63CMMXWw5NVMVdXKknVW07ClOWFYd/p1yuTej/xUVVvpy5hYY8mrmVZs2EZ5ZbWVvFqYJJ9w9eH9KVhXyr++q935gTHGa01KXiJytojME5Ft7t8cEfl9pIKLBzUdUFqbhi3PhL27MKR7G+7/2E95ZbXX4RhjgjQ6eYnI2cAfgWuAbji9IF8PXCkiZ0UkujiQX1SKCPTrZG0atjQ+n3DtEQNYvWk7r3y70utwjDFBmlLyugQ4QVU/VdUtqrpZVT8BTgIui0x4sc9fHKB3h9a0SrWGXVuiQ/rncEBuBx76pIDt5VVeh2OMcTUlebVR1cLaA91hbcIVULxZXLTV7ne1YCLCtUcOoCSwk2dnFXodjjHG1ZTkVV+rpQnZoumOiioKN5RZTcMW7oA+HTikfw7/+HwpW3dY6/PGxIKmJK9BIvJDHX8/AgMjFWAsW1ayjapqtZJXArj2iAFsLqvgyS+Xex2KMYamdUY5KGJRxKldNQ2t5NXi7dOjLUfv04WnvlzG2Qf1pqP122aMp5pS8roe6KGqK+r6i1SAsSy/OEBKkpDbMcPrUEwUXH14f7ZXVDHl86Veh2JMwmtK8vID94hIoYjcJSLDIhVUvPAXBdgrO5PUZHvWOxH065TFCcN68MysFdZxpTEea/RRV1UfVNWDgEOADcDTIrJYRP5XRPpHLMIYll8csEuGCeaP4/NQVR7+pMDrUIxJaE0uMriXCf+uqsOAicDxwKJwBxbrSndWsnrTdkteCaZnh9ZMPKAXr367isL127wOx5iE1eTkJSLJInKMiLwAvAfkAyeGPbIYt8StrGE1DRPP5Yf2IzlJeOBjv9ehGJOwmtI81OEi8jSwGpgEvAP0VdXTVfXtSAUYq6xNw8TVqU0654zqw9vf/0y+24u2MSa6mlLyeh+YCQxS1WNV9UVVTdjrJvlFpbRKSaJH+1Zeh2I8cPEhe5GZmsy9H+Z7HYoxCakpyesHVX1SVTdFLJo44i8O0L9zJj6feB2K8UC71qlMGrMXHy4sZv6qzV6HY0zCacpDyjkicnWokap6XxjiiRv5xQHG9s/xOgzjofMO7sPUmYXc80E+z19woNfhGJNQmlLySgKy6vlLGBu3lVMS2Gk1DRNcZloyl47ty1cF65m5dL3X4RiTUJpS8lqrqn+JWCRxxG81DY3rzJG9eeqr5dzzQT5vXNIREbuMbEw0NKXkZXulq6aGmZW8THpKEn84LI95KzfzyeJ1XodjTMJoSvIaF7Eo4kx+cYC2rVLolGWNsxo4ZUQPendszd0f5FNdrV6HY0xCaErzUBsjGUg88RcFGNA5yy4RGQBSknxcfXh/FhcFeOfHtV6HY0xCiGqLsiIyQUTyRaRARG6oZ7qTRERFZEQ042sMVSW/OED/Lpleh2JiyDH7dmNglyzu+8hPZVW11+EY0+JFLXmJSBIwGTgKGAxMFJHBdUyXBVwJfB2t2JqiaOsOAjsqrWUNsxufT7j68P4sX7+NN+at9jocY1q8aJa8DgAKVHWZqpYDLwPH1THd7cDfgR1RjK3RaiprWE1DU9vhgzsztGc7Hvx4CTsrq7wOx5gWLZrJqzuwKuj9anfYLiKyP9BTVd+pb0YicqGIzBGROSUlJeGPtB7We7IJRUS4/sgB/LxlBy9+vdLrcIxp0WKmF0UR8QH3Adc0NK2qPq6qI1R1RE5OdFu5yC8qpXObNNq1To3qck18GN0vm1F9OzL50wK27az0OhxjWqxoJq81QM+g9z3cYTWygCHAZyJSCIwEpsVapQ2nTUMrdZnQrj1yAOtLy5k6s9DrUIxpsaKZvL4F8kSkj4ikAqcD02pGquoWVc1W1VxVzQVmA8eq6pwoxlivqmplybqAVdYw9dq/V3vGD+rElM+XsqWswutwjGmRopa8VLUSuBz4AKfn5VdVdYGI3CYix0YrjuZYtbGMHRXV9Lf7XaYB1xwxgMCOSh77YqnXoRjTIjWlbcNmU9V3gXdrDbslxLRjoxFTU+RbB5SmkQZ1bcOxQ7vxzxmFnDu6DznWGosxYRUzFTbigd+tJp/X2R5QNg276vD+lFdVM/nTAq9DMabFseTVBPnFAXp1aE3r1KgWWE2c6pOdwSnDe/Di1ytZvanM63CMaVEseTWB1TQ0TXXFuDwAHpq+xONIjGlZLHk1UnllNctKtjHA2jQ0TdCtXSvOHNmbN+atYWlJqdfhGNNiWPJqpGXrS6msVit5mSa79NC+pCX7uP8jv9ehGNNiWPJqJOuA0uyp7Mw0zj+4D//5YS0Lft7idTjGtAiWvBrJXxwg2SfslW2XDU3TXfDrvWiTnsy9H1rpy5hwsOTVSPlFpfTJziA12TaZabq2rVK4eGxfPlm8jrkrrF9XY5rLjsSN5C8OWMsaplnOGZVLdmYad72fj6p6HY4xcc2SVyOUlVeycmOZtaxhmqV1ajJ/OKwfXy/fyFcF670Ox5i4ZsmrEZYUO1WcrbKGaa7TD+hJ93atuPsDK30Z0xyWvBrB2jQ04ZKWnMSV4/P4YfUWPlhQ7HU4xsQtS16N4C8KkJ7io2eH1l6HYlqAE4d1Z6+cDO77KJ+qait9GbMnLHk1Qn5xgLxOWST5xOtQTAuQnOTjmsMH4C8uZdr3axr+gDHmFyx5NYK1aWjC7aghXdi7Wxvu/2gJ5ZXVXodjTNyx5NWAzWXlFG/daW0amrDy+YRrjxjAyo1lvDpnldfhGBN3LHk1wO/WNLSSlwm3sQNyGNG7PQ9/soQdFVVeh2NMXLHk1YBdNQ2tmrwJMxHhuiMHULx1J8/NWuF1OMbEFUteDfAXBchKT6ZLm3SvQzEt0IF7dWRM/xwe/ayAwI4Kr8MxJm5Y8mpAflGAAZ2zELGahiYyrjtiAJvKKnjqq+Veh2JM3LDkVQ9VJd/aNDQRtk+PtkzYuwtPfrmcTdvKvQ7HmLhgyase6wI72bK9wlrWMBF3zRH92VZeyZTPl3odijFxwZJXPWo6oLSahibS8jpnccKw7kydWUjx1h1eh2NMzLPkVQ9/cU3ysme8TORdNb4/1ao8/MkSr0MxJuZZ8qpHflGA7Mw0OmameR2KSQA9O7TmtF/15OVvVrFyQ5nX4RgT0yx51cNfHGCgVdYwUfSHw/JI8gkPTPd7HYoxMc2SVwjV1Yq/uNTud5mo6twmnXNG5fLWd2t2XbY2xvySJa8QVm/azvaKKmvT0ETdxYf0JSM1mfs+tNKXMaFY8gohv9hqGhpvtM9IZdKv9+L9BUX8sHqz1+EYE5MseYVQc8kmz5KX8cB5B+fSvnUK91jpy5g6WfIKIb8oQI/2rchMS/Y6FJOAstJTuHRsP77wl/C39xZRbT0uG7MbOzKH4C8OWMsaxlPnjs5lxcZtPPb5MpaVbOOB0/Yjw06mjAGiXPISkQkiki8iBSJyQx3jrxaRhSLyg4hMF5He0YyvRkVVNUtLSq1NQ+Op5CQftx83hL8cuzfTFxVz8pRZrNm83euwjIkJUUteIpIETAaOAgYDE0VkcK3JvgNGqOq+wOvAXdGKL1jh+m1UVKmVvIznRISzR+Xyz3MPYPXGMo57ZAbzVm7yOixjPBfNktcBQIGqLlPVcuBl4LjgCVT1U1WtaVpgNtAjivHtstjaNDQx5pD+Obx56ShapyZx+uOzeXv+Gq9DMsZT0Uxe3YFVQe9Xu8NCOR94L6IRheAvDpDkE/bKyfBi8cbUKa9zFv+6bDT79WzHlS/P576P/FaRwySsmKxtKCJnAiOAu0OMv1BE5ojInJKSkrAvP78oQG7H1qSnJIV93sY0R4eMVJ4//0BOGd6Dh6Yv4Q8vfcf28iqvwzIm6qKZvNYAPYPe93CH7UZExgM3Aceq6s66ZqSqj6vqCFUdkZOTE/ZA/cUBBlhlDROjUpN93HXyvtx09CDe/Wktpz0+y7pRMQknmsnrWyBPRPqISCpwOjAteAIRGQY8hpO41kUxtl22l1exYmOZ3e8yMU1EmDRmL544awRL15Vy7CNf8dOaLV6HZUzURC15qWolcDnwAbAIeFVVF4jIbSJyrDvZ3UAm8JqIzBeRaSFmFzEF60pRxWoamrgwfnBnXr9kFMk+HydPmcl7P671OiRjoiKqTzyq6rvAu7WG3RL0enw046lLTZuGdtnQxItBXdvwr8tGc9Fzc7jkhXlce0R/Lju0HyLidWjGRExMVtjwkr84QGqyj94draahiR85WWm8OGkkx+/XjXs+9HPVK/PZUWEVOUzLZW3N1JJfFCCvUyZJPjtrNfElPSWJ+0/bj7zOWdz9QT4rN5bx2FkjyMmynsBNy2Mlr1qsTUMTz0SEyw7txz9+tz8L127l+MkzWLR2q9dhGRN2lryCbNlewdotO6xNQxP3jtqnK69fPIrK6mpO/sdMPl5Y7HVIxoSVJa8gS2oqa1jJy7QAQ7q3ZdrlB9O3UyaTnpvD418sRdVa5DAtgyWvILt6T7aSl2khOrdJ55ULD+LoIV25893F/OmNHyivrPY6LGOazSpsBPEXBchMS6Zb23SvQzEmbFqlJvHwxGH07ZTJQ9OXULihjClnDqdDRqrXoRmzx6zkFSS/OED/zpn2fIxpcXw+4erD+/Pg6fsxf9Vmjp88Y9dlcmPikSUvl6qSX2RtGpqW7bj9uvPKhSMpK6/ixEdn8lm+J62wGdNslrxcJaU72VRWYW0amhZvWK/2vH35aHp0aM15U79l6ozlVpHDxB1LXi5/USlgNQ1NYujerhWvX3wQ4wZ15tZ/L+R/3v6JiiqryGHihyUvl9U0NIkmIy2Zx84czsWH9OX52Ss555/fsKWswuuwjGkUS14uf1GAjhmpZGdaUzomcfh8wg1HDeTuk/flm+UbOeHRGSxfv83rsIxpkCUvl1PT0EpdJjGdMqInL04ayebtFRw/eQYzC9Z7HZIx9bLkBVRXK0us92ST4H6V24G3LxtN5zZp/P7pb3jx65Veh2RMSJa8gDWbt7OtvMqSl0l4PTu05o1LRnFwXjZ/futH/vLvBVRaRQ4Tgyx54bQkD9hlQ2OArPQUnjr7V5w3ug//nFHIBc/OYesOq8hhYoslL4JqGnbO9DgSY2JDkk+45ZjB3HnCPny1ZD3HPzKDjxYW2/NgJmZY8sKpadi9XSuy0lO8DsWYmHLGgb147vwDAZj07BxOmTKLbws3ehyVMZa8AMgvLrVSlzEhHNS3Ix9eNYa/nbgPqzaVccqUWZw39Vvr5NJ4KuGTV2VVNUvXldrDycbUIznJx8QDevHZtYfypwkDmVO4kaMf+pKrXpnPqo1lXodnElDCJ6/CDWWUV1Vbs1DGNEKr1CQuGduXL68/jIvG9OXdH9dy2L2fceu0Bawv3el1eCaBJHzyspqGxjRd29Yp3HDUQD6/7lBOHt6T52avYMxdn3LfR34CVjPRREHCJ6/8ogA+gX6d7J6XMU3VpW06fztxHz66agyHDujEQ9OXcMjdn/HUV8vZWVnldXimBbPkVRQgt2MG6SlJXodiTNzaKyeTyb/bn2mXj2ZQ1yxu/89CDrvnc96Yu5qqaqteb8Iv4ZOX39o0NCZs9u3RjhcuGMnz5x9Ih4xUrnnte45+8Es+tmfETJgldPLaUVFF4YZtVtPQmDA7OC+bty8bzeQz9qe8qpoL7BkxE2YJnbwK1pVSrdYBpTGR4PMJv9m3Kx9eNYY7T9iHlRudZ8TOn/oti4vsGTHTPAmdvGpqGg7oYpU1jImUlCQfZxzYi8+vc54R+7ZwI0c9+CVX2zNiphkSOnnlFwdITfLRu2OG16EY0+IFPyN24Zi9eMeeETPNkNDJy18UoG+nTFKSEnozGBNVbVuncONRg9xnxHrw3OwVHHLXp9z/kZ/SnZVeh2fiREIftf3FpQywNg2N8YTzjNi+fHjVGA4ZkMOD05cw5q5PedqeETONkLDJK7CjgjWbt1tNQ2M81jcnk0d/N5y3LxvNwC5Z3GbPiJlGiGryEpEJIpIvIgUickMd49NE5BV3/NcikhupWPzFpYDVNDQmVgzt2Y4XJ43kufMPoH1Gyq5nxKYvsmfEzC9FLXmJSBIwGTgKGAxMFJHBtSY7H9ikqv2A+4G/Ryoea9PQmNj067wcpl12MI+cMYydlVWc/4w9I2Z+KTmKyzoAKFDVZQAi8jJwHLAwaJrjgFvd168Dj4iIaAROu/bp3pY/js+je7tW4Z61MaaZfD7ht/t248i9u/DqnFU8+PESTpkyiz7ZGST7xOvwYsZFh/Tl5OE9vA7DE9FMXt2BVUHvVwMHhppGVStFZAvQEVgfPJGIXAhcCNCrV689CmZI97YM6d52jz5rjImOlCQfvzuwNycO68Ezswr5YfVmr0OKKe1aJW7v79FMXmGjqo8DjwOMGDHCLoYb08K1Sk3i4kP6eh2GiSHRrLCxBugZ9L6HO6zOaUQkGWgLbIhKdMYYY+JGNJPXt0CeiPQRkVTgdGBarWmmAWe7r08GPonE/S5jjDHxLWqXDd17WJcDHwBJwNOqukBEbgPmqOo04CngOREpADbiJDhjjDFmN1G956Wq7wLv1hp2S9DrHcAp0YzJGGNM/EnYFjaMMcbEL0texhhj4o4lL2OMMXHHkpcxxpi4I/FeE11ESoAVe/jxbGq13pHgbHvszrbHf9m22F1L2B69VTXH6yD2VNwnr+YQkTmqOsLrOGKFbY/d2fb4L9sWu7Pt4T27bGiMMSbuWPIyxhgTdxI9eT3udQAxxrbH7mx7/Jdti93Z9vBYQt/zMsYYE58SveRljDEmDlnyMsYYE3cSNnmJyAQRyReRAhG5wet4vCQiPUXkUxFZKCILRORKr2Pymogkich3IvIfr2Pxmoi0E5HXRWSxiCwSkYO8jskrInKVu4/8JCIviUi61zElqoRMXiKSBEwGjgIGAxNFZLC3UXmqErhGVQcDI4HLEnx7AFwJLPI6iBjxIPC+qg4EhpKg20VEugNXACNUdQhO107WbZNHEjJ5AQcABaq6TFXLgZeB4zyOyTOqulZV57mvAzgHp+7eRuUdEekB/AZ40utYvCYibYExOH3toarlqrrZ06C8lQy0cnt6bw387HE8CStRk1d3YFXQ+9Uk8ME6mIjkAsOArz0OxUsPANcD1R7HEQv6ACXAP93LqE+KSIbXQXlBVdcA9wArgbXAFlX90NuoEleiJi9TBxHJBN4A/qiqW72Oxwsi8ltgnarO9TqWGJEM7A/8Q1WHAduAhLxHLCLtca7Q9AG6ARkicqa3USWuRE1ea4CeQe97uMMSloik4CSuF1T1Ta/j8dBo4FgRKcS5nHyYiDzvbUieWg2sVtWakvjrOMksEY0HlqtqiapWAG8CozyOKWElavL6FsgTkT4ikopz03WaxzF5RkQE557GIlW9z+t4vKSqN6pqD1XNxfldfKKqCXt2rapFwCoRGeAOGgcs9DAkL60ERopIa3efGUeCVl6JBcleB+AFVa0UkcuBD3BqDD2tqgs8DstLo4GzgB9FZL477M+q+q53IZkY8gfgBfdEbxlwrsfxeEJVvxaR14F5ODV0v8OaifKMNQ9ljDEm7iTqZUNjjDFxzJKXMcaYuGPJyxhjTNyx5GWMMSbuWPIyxhgTdyx5mbgmIlUiMl9EvheReSISlodGRWRsY1qUF5HPRGSE+/pdEWkXjuU3hYjcJiLjo71cY7yUkM95mRZlu6ruByAiRwJ/Aw7xIhBVPdqj5d7ixXKN8ZKVvExL0gbYBE6rISJyt9vv0o8icpo7fKxbWqrpn+oFt7WEmj7eFovIPODEuhYgIq1E5GW3X6u3gFZB4wpFJFtEct35TBURv7uM8SIyQ0SWiMgB7vQZIvK0iHzjNnp7nDv8HBF5U0Ted6e/yx2e5M6zZp2ucodPFZGT3dfj3Hn96M47LSi2v7il0x9FZGBEvgFjosRKXibetXJbBUkHugKHucNPBPbD6X8qG/hWRL5wxw0D9sbpzmIGMFpE5gBPuJ8vAF4JsbxLgDJVHSQi++K0tlCXfsApwHk4zZGdARwMHAv8GTgeuAmn+anz3MuN34jIx+7n93Pj3Anki8jDQCegu9uXFLUvUbodI04FxqmqX0SedeN9wJ1kvaruLyKXAtcCF4SI3ZiYZyUvE++2q+p+bkeJE4Bn3ZLUwcBLqlqlqsXA58Cv3M98o6qrVbUamA/kAgNxGl1dok6zM6Ea4x1TM05VfwB+CDHdclX90V3GAmC6O98f3eUBHAHc4Cbfz3AScC933HRV3aKqO3DaEuyN0zTTXiLysIhMAGq3/D/AXa7fff+MG2+NmgaX5wbFYExcsuRlWgxVnYVTysppYNKdQa+riMwViOBlVAe9rw5angAnucl3P1XtpaqL6vh8FZCsqptwSpKfARfT9M4ya+YZqXU2JmoseZkWw72PkwRsAL4ETnPvE+XglEC+qefji4FcEenrvp8YYrovcC4BIiJDgH2bEfIHwB+C7rkNq29iEckGfKr6BnAzv+yaJB9nHfq578/CKXEa0+LY2ZeJdzX3vMApyZytqlVuZYqDgO8BBa5X1aJQFRVUdYeIXAi8IyJlOMkvq45J/4HTq/AinO4wmtNp5e0496N+EBEfsBz4bT3Td3eXXXPSeWMd63Au8Jo43dR/C0xpRnzGxCxrVd4YY0zcscuGxhhj4o4lL2OMMXHHkpcxxpi4Y8nLGGNM3LHkZYwxJu5Y8jLGGBN3LHkZY4yJO/8P/D7epgNPTQ0AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "TVD = []\n", "for chi in range(1,11):\n", " mps.set_cutoff(chi)\n", " mps.set_input_state(input_state)\n", " probs_mps = mps.prob_distribution()\n", " TVD.append(tvd(probs_mps, probs_slos))\n", "\n", "plt.plot(TVD)\n", "plt.xlabel('Bond dimension')\n", "plt.ylabel('TVD')\n", "plt.title('Evolution of the TVD between SLOS and MPS probability distributions');" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the TVD decreases as the size of the matrices increases, untill reaching 0 for a bond dimension of 7." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "\n", "> [1] Ulrich Schollwöck. The density-matrix renormalization group in the age of matrix product states. [Annals of Physics](https://doi.org/10.1016/j.aop.2010.09.012), 326(1):96–192, jan 2011.\n", "\n", "> [2] Changhun Oh, Kyungjoo Noh, Bill Fefferman, and Liang Jiang. Classical simulation of lossy boson sampling using matrix product operators. [Physical Review A](https://doi.org/10.1103/PhysRevA.104.022407), 104(2), aug 2021.\n", "\n", "> [3] Hui Wang, et al. Boson Sampling with 20 Input Photons and a 60-Mode Interferometer in a $10^{14}$-Dimensional Hilbert Space. [Physical Review Letters](https://link.aps.org/doi/10.1103/PhysRevLett.123.250503), 123(25):250503, December 2019. Publisher: American Physical Society." ] } ], "metadata": { "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 2 }